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The analytical expression for angular integral flux of atmospheric muons in matter with the 
explicit relation of its parameters with those of the sea level spectrum is obtained. The fitting formula 
for the sea level muon spectrum at different zenith angles for spherical atmosphere is proposed. 
The concrete calculations for pure water are presented. Fluctuations of muon energy losses are 
taken into account by means of parametrized correction factor calculated using survival probabilities 
resulted from Monte Carlo simulations. Parametrizations of all continuous energy losses are obtained 
with using the most recent expressions for muon interaction cross-sections. The corresponding 
parametrization errors and field of method application are comprehensively discussed. The proposed 
formulae could be useful primarily for experimentalists processing data of arrays located deep under 
<^> water or under ice. 
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I. INTRODUCTION 



The last several years have been marked by starting of fullscale data taking of large neutrino and muon telescopes 
located at lake Baikal (NT-200 ||) and in deep polar ice (AMANDA J|). Two underwater telescopes ANTARES || 
' and NESTOR Q] assuming installation at bigger depths are intensively constructed in the Mediterranean. The 
possibility of deployment of telescopes with huge detecting volumes up to 1 km 3 (KM3) is also under wide investigation. 

So, the knowledge of expected angular distribution of integral flux of atmospheric muons deep underwater is of 
interest not only for cosmic ray physics but also for estimation of possible background for neutrino detection and at 
last for a test of the correctness of underwater telescope data interpretation by using the natural flux of atmospheric 
muons as calibration source. The last item frequently implies the estimation with an appropriate accuracy (e.g., 
better than 5 % for a given sea level spectrum) the underwater integral muon flux for various sets of depths, cut-off 
energies and angular bins especially for telescopes of big spacial dimensions (e.g., AMANDA, KM3). 

The existing methods of calculations of muon propagation through thick layers of matter (basically for a standard 
' rock) are largely presented in the literature. There are two completely different approaches in such calculations. In 
first (historically) works on the subject of muon propagation the integro-differential kinetic equation for muon flux 
i-C ■ had been formulated and approximately solved using semianalytical methods. We mention here one of the pioneer 
works of this approach: Ref. |5|, and one the most recent papers: Ref. |(| (this latter paper contains also the large 
bibliography). Works of the second approach use Monte Carlo (MC) technique for propagation studies: the main 
element in such works is simulation of sequences of free flights and interactions during muon passage through a 
medium. One of the most early works using this method is Ref. 0, the most recent one is the work of the present 
authors Ref. ||, containing also many bibliographic references. 

Up to now the presentation of the results of these calculations both for parent muon sea level spectra (especially for 
angular dependence taking into account the sphericity of atmosphere) and for underwater angular flux was not quite 
convenient when applying them to concrete underwater arrays. In addition, a part of numerical results is available 
only in data tables (often insufficient for accurate interpolation) and figures. The possibility of direct implementation 
of MC-methods depends on the availability of corresponding codes and usually assumes rather long computations and 
accurate choice of the grid for simulation parameters to avoid big systematic errors. 

Therefore, the necessity of analytical expressions both for underwater muon integral flux and underwater muon 
differential spectrum is still actual. Besides, the possibility of reconstructing the parameters of a sea level spectrum by 
fitting measured underwater flux in the case of their direct relation looks rather attractive. The only (to our knowledge) 
existing analytical description of muon underwater angular flux published in Ref. (^] is rather comprehensive but its 
relation to widely used parametrizations of sea level spectrum is unclear. 

In this paper we present rather simple method allowing to calculate analytically the angular distribution of integral 
muon flux deep under water for cut off energies (1-10 4 ) GeV and slant depths of (1-16) km for conventional (it, K) sea 
level spectra of atmospheric muons having widely used parametrization by means of 5 parameters. The fluctuations 
of muon energy losses are taken into account. 

The plan of the paper is as follows. In Sec. || the idea of the method used for the derivation of the analytical 



expression for an underwater integral muon flux is described. In Sec. HI the common analytical description of widely 
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used sea level muon angular spectra is introduced. In Sec. [V the parametrization of muon continuous energy losses 
is given. In the same section the calculation of survival probabilities needed for a taking into account muon energy 
loss fluctuations is discussed. In Sec. [v] the concept of an correction factor is introduced. The results of numerical 
calculations and analytical parametrization of this factor are given. At the end of this section the final analytical 
expression for the angular integral underwater flux is derived. The main conclusions of the work are given in Sec. [v]. 
Some mathematical details are presented in Appendices ^ and [b|. 



II. METHOD 



The integral flux of muons in matter with energies above cut-off Ef expected at slant depth R at a zenith angle of 
taking into account fluctuating character of muon losses is conventionally described by 



/>oo 

Ifi(>E f ,R,6)= / P(E ,R,>E f )D(E ,e)dE , 
Jo 



(2.1) 



where P(Eq, R, > Ef) is the probability that a muon, having starting energy Eq, after passing of path R will survive 
with final energy above cut-off Ef and D(Eq,9) = dN(Eo,8)/dEo is a sea level differential angular spectrum. In a 
case of a flat surface the slant depth is expressed by R = h/ cos 6, where ft, is a vertical depth below the surface. 

With the assumption of continuous energy loss rate of muon in matter, L(E) = —dE/dx = a(E) + b(E)E, the 
integral flux I c i is derived from the equation 

/>oo 

Ici(>E f ,R,0)= / D(E ,6)dE , (2.2) 

J E s 

where the low limit of the integral, E Sl is the value of starting energy Eq which results, after passing of path R, in 
the final energy Ef. This value is derived from the solution of the integral equation 

ES dE R. (2.3) 



E, 



L{E) 



The method proposed in this work for an obtaining the approximate analytical expression for the real integral muon 
flux allowing for loss fluctuations consists in the following. Firstly, we calculate numerically Ifi and I c i and the 
correction factor defined by the ratio 

In principle, this factor can be calculated using known codes for muon propagation through matter. In this work we 
apply for this aim the MUM code described in previous work of the present authors |§. Further, we approximate the 
muon energy loss function L(E) by a linear energy dependence, 

L(E)=a + f3E (2.5) 

(a and are independent on energy) and parametrize the sea level muon spectrum by dependencies of the type 

D(E ) oc -'< 



1 + E /e Q ' 

Using these approximations, the integral flux in matter defined by Eq. ( |2.2| ) can be expressed analytically. Below in 
this paper we will designate this integral flux by F c i (to distinguish it from I c i calculated numerically). 

The real integral flux allowing for energy loss fluctuations can now be obtained using the correction factor introduced 



in Eq. (2.4). As it is shown below, the dependencies of this factor on its arguments are rather smooth and can be easily 



parametrized, so the final result for the real integral muon flux in a matter is also given by the analytical expression: 
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III. SEA LEVEL SPECTRA OF ATMOSPHERIC MUONS 



From the large collection of conventional (tt, K) ang ular spectr a of atmospheric muons at sea level D(Eq, 9) proposed 



in literature |10|-|19|], for the numerical calculations of (|2.l|) and (2.2) we have chosen as basic the spectrum from work 



of Ref. which had been obtained using calculations described in Refs. JUTLM. The reasons are the following: 

(i) the energy range of its validity is up to 10 9 GeV allowing the calculation of underwater flux up to cut-off energy 
greater than 10 TeV for slant depths as high as (20-30) km; 

(ii) access to data tables (given to us by authors of work [jll[) detailed enough to perform splines via energy and 
zenith angle with acuracy < 0.5%; 

(iii) high level of the agreement of this spectrum with full set of existing experimental results. 

This spectrum results from the computations based on nuclear cascade model of Refs. |2l],^(| (see also Ref jnj) and 
semiempirical model for primary spectrum proposed in Ref. [ p2| . 

Following Ref. |2j| , for the description of basic spectrum for the energy range Eq < 1 .53 TeV we have used values defined 
by D(E ,9) = D(E ,0°)A(E ,9) with vertical spectrum D(E o ,0°) taken from Ref. |l^] and angular distribution 
A(Eq,6) = D(Eo,0)/D(Eq,O°) from data tables based on calculations of Ref. jfl). For energies above 1.53 TeV we 
have used values of D(Eq, 9) according to splines computed via the same data tables. For the spectrum composed by 
this way the abbreviation "NSS" (following names of authors of Ref. fl2(j]) is used below. 



The most convenient analytical form of D(Eo, 9) has been proposed in Refs. 24 25J] and became almost conventional 
after work |j~7|]. It is as follows: 

dN(Eo,8) „ , / 1 By 



dE - D{Ea > 0) ~ Aq {i + EtfE~ (9) + l + Eo/E- (9) ) ^ (3 ' 1} 

Here, 7 is a spectral index and the functions Eq t (9) are given by the relation 

E c : K (9)=E c : K (0°)/cos8*. (3.2) 

These functions have approximate sense of critical energies of pions and kaons for given zenith angle and Eq~ ^(0°) 
are those for vertical direction. The value cos 9* has sense of effective cosine taking into account the deviation of real 
atmosphere from the flat one. Thus, sea level spectrum may be expressed by 5 parameters (A , B , E§ r (0°), Eq' (0°), 7) 
each of them having the physical meaning. 



As it is shown in Appendix H, just parametrization (3.1) is very convenient for a derivation of an analytical 



expression for integral flux in matter. To check the validity of this analytical expression for the large range o f ze nith 



angles and cut-off energies we have performed the fit of original NSS spectrum by 5 parameter expression (3.1) for 
vertical direction. This fit was done within ranges Eq= (0.2-200) TeV. 

When checking the values of fit spectrum for cos 9=(Q. 05-1.0) we realized that the standard description of effective 
cosine (with geometry of spherical atmosphere and with definite value of effective height of muon generation) is not 
enough and one should introduce an additional correction S(9) leading to (10-20) % increase of effective cosine value 
for cos 9 < 0.1. The reason of an appearing of this correction is that the concept of an effective generation height is 
approximate one. It fails at large zenith angles where the real geometrical size of the generation region becomes very 
large. 

The resulting fit of angular sea level spectrum in units of (cm _2 s _1 sr _1 Gey _1 ) is given by 



D(E ,9) = 0.175 
with modified effective cosine expressed by 



0.037 



V 1 



£V)C0S6>** En cos 



q** I ^0 



Eq • , (3.3) 



1 



103 810 



cos0** = S(6) cos 6*, (3.4) 



where cos 9* is derived from spherical atmosphere geometry and is given by formulae of Appendix ^| and S(9) is the 
correction which is given for sec 9 < 20 by 

S(9) = 0.986 + 0.014 sec 8. (3.5) 
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The dependencies of both modified and geometrical effective cosines as function of zenith angle are presented in Fig. [|. 
Note that for cos 6* > 0.4 the influence of the curvature of real atmosphere is less than 4 % but for cos 9 < 0.1 it is 
greater than 40 %. 



We should note that our expression (3J) for the sea level muon spectrum does not contain a contrubution from 
atmospheric prompt muons. According to the most recent calculations based on perturbative QCD, this contribution 
becomes essential only at Eq > 10 6 GeV p6| . The predictions of nonperturbative models (see, e.g., ^7j) are slightly 
more optimistic. We plan to generalize our approach and include this contribution in our following paper. Incidentally, 
inclusion of prompt muons should be done in parallel with taking into account the steepening of the sea level muon 
spectrum due to the knee in the primary cosmic ray spectrum. 



Fig. H illustrates the limits of applicability of angular spectrum given by Eq. (3.3), for energy and zenith angle 



variables. The energy region, inside which the deviation from parent NSS spectrum is less than 5 %, is shifted from 



(0.3-200) TeV for cos (9=1.0 to (1.5-300) TeV for cos 0=0.05. The sea level spectrum given by (^3|) is valid only below 
the knee (Eq ~ 300 TeV) of primary cosmic ray spectrum. Below we will check the influence of the knee to the value 
of underwater flux by comparison with calculations resulted from NSS sea level spectrum which takes into account 
the knee. 

Despite the parametrization ( |3.3| ) was done using absolute values of D(Eq,6), it reveals rather small deviation of 
the angular distribution A(E ,8) = D(E ,0)/D(E ,0 o ) (shown in Fig. ||) from that of parent NSS spectrum. Note 
rather good consistence of angular distribution of NSS spectrum with that of Ref. jl5j for cos 9 < 0.2 and differences 
as big as 10 % from that of Ref. |2j| for energies E > 10 TeV and cos# < 0.2. The angular distributions for E > 100 
TeV are not shown in Fig. because of their coincidence with those for £?q =10 TeV. 



The widely used vertical sea level spectra which can be parametrized by (3.1) are presented in Table ||. For further 



TABLE I. Parameters of the fitting formula ( fc.l| ) for some selected vertical sea level muon spectra used in this paper. 
The values of critical energies are in (GeV) and values of other parameters allow the calculation of spectrum in units of 
(cm~ 2 s _1 sr _1 Gel /_1 ). Corresponding energy ranges of fit validity and references are mentioned. 



Name 

VZK 

Gaisser 

Present work 

LVD 

MACRO 



_Ao 

0.1258 

0.14 

0.175 

0.2576 

0.26 



_Bo 

0.0588 

0.054 

0.037 

0.054 

0.054 



#£1(0°) 

100 

104.545 
103 

104.545 
104.545 



Eo r K (0°) 
650 

772.727 
810 

772.727 
772.727 



2_ 

2.65 
2.70 
2.72 
2.77 
2.78 



Eq (TeV) 
0.1 - 100 
0.03 - 100 
0.3 - 250 
2.0 - 40 
0.5 - 20 



Ref 

|1E| 

\m ,_, 

Eg- (ty) 
\m 

I IS I 



investigations we have used most rec ent L VD spectrum presented in Ref. jl8| . The original spectrum of Ref. |l6| (VZK) 
can not be directly decomposed to (3.1) because A$ depends logarithmically on energy. Thus, we were managed to 
construct a "simplified" VZK spectrum with Aq = 0.1258 fc, where k is the coefficient depending on the slant depth 
R, for which underwater flux is calculated, approximately, as: A:=0.87 for R=(l-7) km, k—0.79 for i?=(7-ll) km and 
fc=0.72 for R > 11 km. 



The expression (3T) can be rewritten in the shorter form which will be needed for later sections: 



D(E ,( 



B7-7 
^0 



E 



D 0i 



i + E /E%(ey 



(3.6) 



where Z?o x = Aq and Dq k = AqBq. 



IV. MUON PROPAGATION THROUGH WATER 



A. Continuous energy losses 

For the description of continuous energy losses of muon in water needed for solving the integral equation ( [2.3] ) we 
have done parametrizations based on output data from MUM code. All formulae for cross sections used in MUM 
code for muon energy loss computations are described in details in Ref. Q . 

The energy loss rate per unit of path x is given, conventionally, by the expression L(E) = —dE/dx — a(E) + b(E)E, 
where a(E) is loss due to ionization and b(E) = b p (E) + bt,(E) + b n (E) is the sum of coefficients for all radiative 
processes: e + e~ pair production (b p ), bremsstrahlung (bb) and photonuclear interaction (&„). 
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Ionization loss was taken as composition of 2 processes (Fig. a{E) — a c (E) + a e (E), where a c is classic ioniza- 
tion calculated using Bethe-Bloch formula (Refs. |^8|, 29| )and a e results from e-diagrams for bremsstrahlung being 
treated as a part of ionization process following Refs. [3^,0 (7-quantum is emitted by atomic electron). Taking into 
account the last process leads to 1.8%, 3.4%, 5.5% increase of ionization loss for 100 GeV, 1 TeV, 10 TeV muons, 
correspondingly. Applying designations used in Ref. p2| the approximation formula for a c is given by 

a c {E) = a C0 + a ct In ( j W max = — ^— — , (4.1) 

where W max is maximum energy transferable to the electron and , m e are the rest masses of muon and electron. 
The set of cofficients, in units of (10 _3 GeVcm 2 g _1 ), 

a Co = 2.106, a Cl = 0.0950 for E < 45 GeV, 
a C0 = 2.163, a Cl = 0.0853 for E > 45 GeV, 



gives the error for parametrization ( J4. l| ) smaller than 0.2 % for (1-10 ) GeV range. For a e the following polynomial 
approximation [E in units of (GeV)): 

a e (E) = 3.54 + 3.7851n£+1.151n 2 £ + 0.06151n 3 £; (lO^GeVcmV 1 ), (4.2) 



has the error < 0.2 % for (50-10 8 ) GeV range. Finally, the sum of (4.1 ) and (|4.2|) has the error of a(E) approximation 
smaller than 0.2 % for (1-10 8 ) GeV range. 

Radiative energy losses of muon in water has been calculated using the cross sections from the following works: 
Andreev, Bezrukov and Bugaev (Refs. |3^,|33|]) for bremsstrahlung, Kokoulin and Petrukhin (Refs. ^,Q) for direct 
e + e~ pair production, Bezrukov and Bugaev (Ref. |3q] ) for photonuclear interaction. Fig. || sketches corresponding 
b-terms of radiative energy losses. 

In principle MUM code gives values for all b-terms up to 10 9 GeV. Note that the logarithmic rise with an energy 
of the photoabsorption cross section, used in the model of Ref. p5| and resulting in corresponding increase of the 
6„-term, is experimentally proved already up to the photon energy ~ 10 5 GeV. 

For parametrization of b-terms we have used the same functional form as in Ref. |3^| but increased the power 
of polynomial to improve accuracy and enlarge the range of application up to 100 PeV without division on energy 
subintervals (E in units of (GeV)): 
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bi(E) = V" b i:j In 3 E, where i = p, b, n. (4.3) 



3=0 



Corresponding coefficients of this decomposition are collected in Tab. |l| This fit works with typical errors ±(0.2- 
0.4)% within (50-7xl0 7 ) GeV for all b-terms. 
The sum b(E) of fits for b-terms is valid for the energy range (50-10 8 ) GeV with the relative error < 0.5 %. 



TABLE II. Coefficients b tj (10~ 6 cm 2 g _1 ) of the fitting formula (4.3) for b-terms of muon energy losses in water. Maximum 
absolute values of relative errors are shown. 



b-term 


subscript i 


bio 


bn 


U2 


bi3 


bi4 


Max.err,% 


e + e~ pairs, b v 


P 


-11.31 x 10" 1 


7.876 x 10 _i 


-8.192 x KT 2 


3.763 x 10- 3 


-6.437 x 10" 5 


0.8 


bremsstrahlung, bb 


b 


-1.149 x 10" 1 


2.963 x 1CT 1 


-2.165 x 10" 2 


5.630 x 10~ 4 


-2.119 x 10" 6 


0.9 


photonuclear, b n 


11 


3.903 x 1CT 1 


9.355 x Id -3 


-3.378 x 1(T 3 


4.913 x 10~ 4 


-1.216 x 10" 5 


0.6 



Fig. shows muon energy loss rate in water. Note that the radiative losses account for approximately 5 % of 
the value of total losses at muon ener gy 70 GeV, 50 % at 1 TeV and dominate as 95 % at 20 TeV. The sum of all 
parametrizations (|4.l|), (fb^) and (4J3) results in the description of total energy losses L(E) with excellent accuracy 
< 0.3% with varying sign of error for energy range (1-10 8 ) GeV mainly because of b-term error compensations 
(even for region (1-50) GeV). Thanks to th is fa ct jus t th is parametrization of L(E) has been used for the numerical 
computation of integral flux I c i defined by ( |2.2| ) and (2.3). 

Further, we realized that the total losses L(E) may be described with an accuracy better than 2.5% (with varying 
sign of error) for energy range (3-3 xlO 6 ) GeV by 3-slope linear fit 



L(E) =a + (3E, 



(4.4) 
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where 



a = a Q = 2.30, f3 = fa = 15.50 for E < E m = 30.0 GeV, 

a = a x = 2.67, (3 = fa = 3.40 for E 01 < E < E 12 = 35.3 TeV, 
a = a 2 = -6.50, (3 = fa = 3.66 for E > E 12 = 35.3 TeV, 

and a are in units of (10 _3 GeVcm 2 g~ 1 ) and (3 in (10 _6 cm 2 g~ 1 ). The energy losses expressed by this parametrizations 
have sense of effective ones, for example in energy region 30 GeV-35 TeV the values a = 2.67 and (3 = 3.40 repre- 
sent effective ene rgy losses due to ionization and radiative processes, correspondingly. The availability of the linea r 



dependence (4.4) leads to the possibility to derive an analytical formula for underwater integral flux (see Sec.[ 

For simplicity, we will use in the following only 2-slope approximation (with (a±,fa) and (0:2, fa), see the end of 
Appendix |b|) . 

B. Survival probabilities 

Probabilities P(Eq, R,> Ef), that a muon having starting energy Eq, after getting over the path R will survive 
with final energy above cut-off Et (so-called survival probabilities) have been calculated by means of conventional 
ratio N(Eq, R, > Ef)/N to t(Eo) m direct MC-simulations by using MUM code with the following inputs. 

(i) In addition to fluctuations of all radiative processes (by using the differential cross sections chosen as described 
in preceding subsection) the knock-on electron production was also treated as stochastical (see Ref. [pp|). The 
interactions with fraction of energy lost less than v CM t =10 -3 were treated as continuous. 

(ii) We used Ntot(Eo)=10 5 simulations for each given value Eq of monoenergetic beam resulting in accuracy better 
than 1% for P(Eq,R, > Ef) >0.1. The corresponding MUM-subroutine which watches the muon energy Ef 
in different points along the path 7? in a single simulation act was involved to accelerate considerably the 
calculations. 



(iii) The grid of simulation parameters (Eq, R, Ef) was sampled to compute underwater angular flux defined by (2.1) 
numerically as: 



I fl (>E f ,R,6)= f " P(E ,R 1 >Ef)D(E ,d)dE (4.5) 

with an systematic error less than 0.5 % at the depth of Baikal telescope location h =1.15 km for both Ef=(l- 
10 4 )GeV and cos0=(O.O495-l.O). It corresponds to slant depths R = h/ cos8 of (1.15-23.23) km. The value 
Eq 1 " 1 of minimum detectable energy was found to be 302 GeV and the value E™ ax needed for avoiding under- 
estimation of integral flux Ifi was taken as large as 36.3 PeV. Thus, the resulted matrix of survival probabilities 
has dimension of (509 x 97 x 81) and is computed for a grid of parameters (Eq, h sec 0, Ef) each of them uniformly 
distributed on: £ =(3.02xl0 2 -3.63x 10 7 ) GeV with step Alog 10 E o =0.01, sec0=(l-2O.2) with step A sec 0=0.2, 
and -E/=(l-10 4 ) GeV with step A log 10 Ef =0.05. Although the chosen grid was optimized for the vertical depth 
of h =1.15 km, the resulted survival probabilities may be used for calculations at deeper vertical depths up to 
slant depths R of 23.23 km. 

The survival probabilities resulted from MUM code are shown in Fig. [?|. The sampled comparison of our results 
with those obtained with widely used codes PROPMU |}7| and MUSIC @ is presented in Fig. §. Note, that 
the greater survival probabilities given by PROPMU for Eq= (300-800) GeV lead to (15-20)% overestimation of 
Ifi(— 10 GeV,R,8) for the slant depths of i?=(l-2)km that is critical for telescopes located at shallow depths (e.g., 
at lake Baikal). The discrepancy with results given by MUSIC code both for survival probabilities and ///(> Ef,R,9) 
is less than 3 % for a wide range of Ef and R. 

Fig. |^ sketches the values of starting muon energy Eq at a sea level which result in passing till underwater vertical 
depth of 1.15 km at various zenith angles and cut-off energies Ef. 

V. UNDERWATER SPECTRA OF ATMOSPHERIC MUONS 
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A. Correction factor 



The influence of fluctuations of muon losses in matter (mainly due to radiative processes) results in that the real 
integral flux Ifi is generally greater than I c i calculated in the approximation of continuous losses. In this work we 
propose to allow for the influence of energy loss fluctuations on the value of angular flux in matter by means of 
correction factor expressed by the ratio 



C f (>E f ,R,t 



Id(>E f ,R, 9) 
I fl (>E f ,R,6) 



(5.1) 



Here Ifi is given by (2.1) and I c [ is given by (2.2). In the assum ptio n when the same differential cross sections are 
used for computing both the numerator and denominator of ratio (5.1), one may expect that this factor depends only 
weakly on sea level angular spectrum. For numerical calculations of the correction factor we used total continuous 
energy losses defined by sum of (4.1), (|4.2D and (4.3), and survival probabilities described in preceding section. 



The value of underwater integral flux for a given slant depth R at zenith angle 9 is greater than the flux calculated 
for the same R for vertical direction as is shown in Fig. It results from angular distribution of sea level spectrum 
shown in Fig. ||. At the same time the values of correction factors calculated for the same slant depth R at vertical 
direction and at zenith angle 9 differ weakly. It is illustrated in Fig. [ll], where one can see that C/(> Ef, R, 0°) differs 
from Cf(> Ef,R,&rccosh/R) maximum on 3.3% for Ef >10 GeV at vertical depth h of 1.15 km. It appears that 
with acceptable accuracy the correction factor depends on slant depth R only, rather than on R and 9 separately. 

The dependencies of correction factor on Ef and R, calculated for sea level spectrum (3.3), represent the set of 
rather smooth curves (shown in Fig. |Tl"|) and it is possible to approximate this factor with formula 



C f (>E f ,R,0)=£(E, 



v logj 



Ef)R l 



(5.2) 



i=0 3=0 



Here cut-off energy Ef is expressed in (GeV) and slant depth R is in (km) with the coefficients c,*j collected in 



Table HI. When using ( p.2|) for cut-off energies Ef <10 GeV one should substitute value of Ef=1 G eV. 

Formula (5.2) can be applied for any geometrical shape of the surface. Right hand side of ( |5.2| ) depends on 9 
because, generally, R = R(9). So, in the particular case of a flat surface the angular dependence of the correction 
factor appears, in our approximation, only through the relation 



R 



cos( 



where h is a vertical depth. 

The accuracy of formula ( p^ ) for _E/ = (1-100) GeV is better than ±2% for slant depths R as large as 22 km and 
is not worse than ±3 % for Ef=l TeV up to R=17 km and for E f =10 TcV up to R=15 km. Fig. [n] shows that for 



TABLE III. Coefficients c, 
trum (p.3j) used in this work. 



of the fitting formula (p.2|) for correction factor calculated for vertical basic sea level spec- 



subscript i qo Cii ca ca eg 

6.3045 x 10 _1 6.6658 x KT 1 -4.5138 x KT 1 1.2441 x KT 1 -1.1904 x 10 -2 

1 2.0152 x lO^ 1 -4.2990 x 10" 1 3.2532 x 10" 1 -1.0265 x 10" 1 1.0751 x 10~ 2 

2 -3.3419 x 10~ 2 5.1833 x 10~ 2 -3.9229 x 10" 2 1.2360 x 10~ 2 -1.2911 x 10~ 3 

3 1.6365 x 10~ 3 -2.3645 x 10~ 3 1.7775 x 10" 3 -5.5495 x 10~ 4 5.7557 x 10~ 5 

4 -2.6630 x 10" 5 3.7770 x 10" 5 -2.8207 x 10" 5 8.7275 x 10" 6 -8.9919 x 10" 7 



Ef < 100 GeV the total energy loss may be treated as quasi-continuous (at level of Cf > 0.9) only for slant depths 
R < 2.5 km but for Ef— 10 TeV the fluctuations should be taken into account at level of 15 % already for slant depth 
as small as R=l km. For slant depths larger than 10 km Ifi and I c i differ more than on a factor of 2. 

The dependence of correcton factor Cf on different sea- level vertical spectra from the Table | is illustrated by 
Fig. |l2|. The correction factors calculated for Ef=10 GeV using sea level spectrum (3.3) with spectral index 7 of 
2.5 and 3.0 (instead of 2.72) differ more than on a factor of 2 starting from slant depth of R—12 km. Nevertheless, 
the values of Cf calculated using sea level spectra having 7=(2.65-2.78) are already within ±5 % corridor. For Ef=l 
TcV this corridor is larger on 2%. This fact results in the possibility to extrapolate the parametrization ( |5.2| ) based 
on sea level spectrum having 7=2.72 to other spectra from Table | at least up to slant depths of (12-13) km without 
introduction of additional spectral corrections. 
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B. Comparison of numerical and analytical calculations 



Taking into account the results of Appendix [b] and preceding subsection we present finally the analytical expressions 
for calculations of underwater angular flux above cut-off energy Ef for a slant depth R = h/ cos 9 seen at vertical 
depth h at zenith angle 8 and allowing for the fluctuations of energy loss: 



F fl (>E f ,R,t 



Fd(>E f ,R, 9) 

c f (>E f ,R,ey 



(5.3) 



where correction factor is expressed by (5.2) and angular flux F c i(> Ef,R,9) based on effective linear continuous 
energy losses a + (3E having 2 slopes, is calculated by the following rule: 



F c i(>E f ,R, 



F cl (>E f ,R,9; ai.ft) for R < R 12 , 
F cl (>E l2 ,(R-R l2 ),6; a 2 ,p 2 ) for R> R 12 . 



(5.4) 



Here Eyi is the energy in the point of slope change from (an, Pi) to (ot 2 , P2) and R\ 2 is the muon path from the energy 
E12 till Ef which is given by 

1 fai+Eiafa 

R12 = — m 

Pi yax+Effa 

As is shown in Appendix [B| the formula for integral muon angular flux in the assumption of linear continuous energy 
losses is as follows: 



F cl (>E f ,R,6; a, 13) 



-PRrf 



£ D Qi E^(6)(E f + y,r 7 (l - Zi) 1-7 8(^,7), 



(5.5) 



,K 



where subscript i stands over both pion (tt) and kaon (K) terms and 

Vi = ~ (1 - e- fiR ) + KM e-? R , * = ^m±!l , 
P E f + yi 



S(z^) = l + —- + 



2z 2 



KM 



K(°°) 

cos 9* 



7 + 1 (7+l)(7 + 2) ( 7 + l)(7 + 2)( 7 + 3) 

The 5 param eter s (Dq^ , D 0k , E^ (0°), Eq~ k (0°), 7) are those of the sea level spectrum ( |Q| ) taking into account the 
notations in (|3.6| ): 

-Do x = A , D 0k = A B . 

The corresponding angular distrubution should be introduced using an analytical description of effective cosine 
cos 9* taking into account the sphericity of atmosphere, ft should be noted that the description of underwater angular 
flux with the 5 parameters of a sea level spectrum gives the possibility of their direct best fit by using the experimental 
underwater distribut ion. 

The flux value in ( |5.5|) is expressed in units of (cm^ 2 s _1 sr _1 ) and all energies are in (GeV), slant depth R in units 
of (gem -2 ), loss terms a and /3 are in unit s of (10 _3 GeVcm 2 g _1 ) and (10 _6 cm 2 g _1 ), correspondingly. Note that 



when calculating correction factor with (5.2) range R is presented in units of ( km 



To simplify the formula we have used only two slopes from approximation (4.4) when substituting in (5.4): (ai=2.67, 
/3i— 3.40 ) an d (a2—— 6.5, /?2=3.66) with i?i2=35.3 TeV. To examine the angular behaviour of a flux given by the 
formula (|5.3|) b y means of the co mpa rison with numerical calculations based on NSS spectrum, we have used (according 
to the fit Q3.3D ), instead of Eq. (3.2), the modified expression for the critical energies E^ (ff) — E^ jc (O°)/cos0**, 
where cos 9** is given by the expression (3.4). The dependence on the sea level spectrum choice was investigated for 
vertical direction only, by u sing the sea level spectra collected in Table |. 

The examination of (5.5) showed rather quick convergence of series S(z,7) with increase of R and Ef. Therefore, 
for the accuracy of F c i computation better than 0.1 % it is quite enough to take only four first terms of this series 
(up to z 3 ) for all values R > 1 km and Ef in (1-10 4 ) GeV. Even using the two terms leads to the accuracy of 1.3 % 
for (R=1.15 km, E f =l GeV) and <0.5% for (R > 2.5 km, E f > 1 GeV). 

Fig. [l^ shows the comparison of underwater angular integral fluxes allowing for loss fl uctu ations at different basic 
depths h (of location of existing and planned telescopes) calculated both numerically (4.5) for NSS spectrum and 



8 



analytically (5.3) for the spectrum given by ( |3.3| ). Note the coincidence of our numerical result obtained using Monte 
Carlo survival probabilities with that calculated using analytical method in the first paper of Ref. Jll| at level of ±5 % 
for basic depth h of 1.15 and 1.61 km for £'/=(10-20) GeV up to cos# =(0.1-0.15). Fig. |lj gives the comparison of 
numerical results and that given by the formula for different vertical sea level spectra for 2 selected depths of interest: 
10 km and 15 km. 

We realized that the error given by formula ( |5.3| ) for all mentioned sea level spectra is within the corridor of ±(4- 
6)% for all cut-off energies i?/ = (l-10 3 ) GeV and slant depths i?=(l-16)km (corresponding angle is expressed by 
cos 6 = h/R for a given vertical depth h). This is proved for h in a range (1-3) km. For bigger c ut-o ffs of Et=(l- 
10) TeV the corridor of errors is ±(5-7) % for i?=(l-13) km. Note that for the sea level spectrum (3.3), just used for 
Cf parametrization, the erro rs a re smaller on 2%. 

Note that the expression (5.3) may be directly used for an ice after substitution R — > R/p, with p being th e ic e 
density and with an additional error of ~ 2 % for a sea water. In spite of seeming complexity of the formulae ( [5.3] ) . 
(5.4) and (5.5) they may be easily programmed. 

The validity of proposed formula up to cut-off energies 10 TeV allows a ca lcul ation of underwater angular differential 
spectrum D(Ef,R,9) by means of numerical differentiation of expression (5.3). It leads to rather appropriate results 
up to slant depths (11-12) km. We illustrate this in Fig. [l5| by comparison the underwater angular spectra calculated 
by numerical differenti atio n of integral fluxes given both by (4J3) for NSS sea level spectrum and analytically (|5.3| ) 
for spectrum given by (3.3). For a vertical depth h—1.15 km it yields to errors ±4 % for i?/=(20-8 x 10 3 ) GeV for the 
angles corresponding to slant depths R = h/ cos (9 of (1-3) km and ±(6-8) % for Ef — (30-5 x 10 3 ) GeV for the slant 
depths (3-12) km. Even for i?=23.2 km the result is still valid within ±10% but for the very narrow energy region 
E f = (90-300) GeV. 



VI. CONCLUSIONS 



The analytical expression presented in this work allows to estimate for fluctuating losses the integral flux of atmo- 
spheric muons in pure water expected for different zenith angles, cos#=(0. 05-1.0), at various vertical depths at least 
of h=(l-3) km for different parametrizations of the sea level muon spectra. The errors of this expression are estimated 
to be smaller than ±(4-6) % for cut-off energies ranged in £'/ = (l-10 3 ) GeV and slant depths in h/ cos #=(1-16) km. 
The main advantage of the presented formula consists in the possibility of the direct best fit of 5 parameters of parent 
sea level spectrum using angular distribution of underwater integral flux measured experimentally at a given vertical 
depth. 

The validity of this analytical expression with an accuracy of ±(5-7) % for £?/=(10 3 -10 4 ) GeV and slant depths of 
(1-12) km gives also the possibility of estimation the angular underwater differential spectrum (by means of numerical 
differentiation) with error smaller than ±(6-8) % for energies of (30-5xl0 3 ) GeV. 

The accuracy of the presented parametrization of the correction factor as a function of Ef and slant depth R is 
rather high and is about ±5 % for all angles and kinds of the sea level spectrum (assuming that the spectral index 7 
is approximately within (2.65-2.78)). It results in the possibility to use it for an estimating numerically from various 
sea level spectra the value of an angular integral flux allowing for fluctuations of losses without direct Monte Carlo 
simulations. 

The results of this work may be used directly to the estimations in ice (substituting a slant depth in ice in units 
of water equivalent) and with additional error ~ 2 % for a sea water. The proposed method may be adapted to 
estimations in rock after corresponding description of the correction factor and continuous effective losses. 
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APPENDIX A: GEOMETRICAL EFFECTIVE COSINE DUE TO SPHERICAL ATMOSPHERE 



As was shown in Ref [B9[ for the spherical isotermal atmosphere the effective cosine is defined by: 



C0Sr " exp(^)(i _ exf(O) V ARe + H) ' (A1) 
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where 



e = 



1 / cos 2 9R E H 



1 



H/ R E \ 



2Hn 



Hn 



2R E . 



(A2) 



Re is the Earth radius (we used mean value 6367.554 km), Hq = 6.4385 km for T= 220° K, erf is the standard integral 
of probability. H has sense of the effective height of muon generation in spherical atmosphere, above and below that 
approximately the same fluxes of muons are generated. For 9 < 70° we have used value H=15 km and for big angles 
(70° — 90°) the approximation obtained in Ref. |39] by using data tables of Ref. Jf0|: 



H = 34 - 10.5 cos 9 - 120 cos 2 9 + 250 cos 3 



(km). 



(A3) 



Since for £ > 3.9 the value of (1 — erf(£)) is equal to with a good accuracy one should develop the expression 
exp(£ 2 )(l — erf(£)) as a series in £ and rewrite (Al) as: 



cos 9* = £ 



2Hn 



E 



(2n-l)!! 



(A4) 



We have used (A4) for cos# >0.1 (that corresponds to £ > 3.1429) and (Al) for cos #=(0-0.1). In spite of rather 
slow convergence of series in (A4), seven its terms are enough for accuracy of 0.1%. For cos 9 > 0.8 the atmosphere 
can be considered as flat one (the deviation is less t hen 0.3%). For an estimation avoiding the special function and 
series calculations we propose the polynomial fit of (Al): 



cos 9* 



E 



c, COS 



(A5) 



with the coefficients of the decomposition assembled in Table IV. The accuracy of (|A5p is much better than 0.3% 



TABLE IV. Coefficients a of the fitting formula (A.5) for effective cosine with the maximum relative errors. 



cos 6 


CO 


Cl 


C2 


C3 


c 4 


Max.err,% 


0=0.002 


0.11137 














0.004 


0.002^0.2 


0.11148 


-0.03427 


5.2053 


-14.197 


16.138 


0.3 


0.2^0.8 


0.06714 


0.71578 


0.42377 


-0.19634 


-0.021145 


0.7 



except the region cos #=(0.3-0.38) where it may reach the value of 0.7%. 



APPENDIX B: DERIVATION OF FORMULA FOR INTEGRAL SPECTRUM IN MATTER 



The solution of (2^3) in the assumption of continuous linear energy losses L(E) = a + f3E is as follows: 



E s (E f ,R) = -((a+0E f y R -a). 



(Bl) 



Here E s is that value of the starting energy Eq of the muon at sea level which results to passing of path R in matter 
with the final energy Ef. Thus the expression (2.2) for the integral flux in matter (after change of designation 
I c i — > F c i) is transformed to 



/>CO 

F c i(>E f ,R,9)= / D(E s (E Q ,R),9)e /3R dE . 



For simplicity let us keep only one term of series (3.6) for the sea level differential spectrum, 

£§ r (0°) 



D(E , 



l + E /E- r (9) : 



E c r (9) 



cos 9* 



(B2) 



(B3) 
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After substitution of (B3) in (B2) one can show that the resulting expression is reduced to 

F cl (> E f ,R,0) = D E^ r {6)e-^ / (E + x)^(E + y)- 1 dE Q , 

Jej 



(B4) 



where x = aCl- e~P R )/(3 and y = x + E§ r (9)e-' 3R . 

Integral in (B4) has the tabulated solution (Ref. El) resulting in the following analytical expression for integral flux 
in matter: 

F cl (> Ef,R, 9) = DoE^^e-^iEf + y)^B( 7 , 1) 2 Fj ( 7 , 7,7+1, z), (B5) 
where B( 7 , 1) is Bessel function and 2F1 (7, 7, 7+1, z) is hypergeometric function of variable z defined by the expression 

z = E?(e)e-P R /(E f +y). 

All the conditions of convergence of (pa) are fulfiled (let us omit here the proof because of long calculations). The 
expressions for special functions in (|B5|) can be reduced to 



where 



B( 7 ,l) = -, 

7 



S(z,7) = 2 Fi(l,l,7+l ) z) = 1 



2 F 1 ( 7 ,7,7+l,z) = (l-z) 1 ^S(2, 7 ), 



2z 2 



6z 3 



7+1 ( 7 +l)(7 + 2) ( 7 +l)(7 + 2)(7 + 3) 

00 / n 

= l + ^n!z n (n(7+i) 

n=l S=i 



(B6) 



Finally, taking into account both (7r,K)-terms of the sea level spectrum (3.6) the expression for the integral angular 
flux in matter in the assumption of continuous losses is as follows: 



F d (>E f ,R,6) 



-mi 



]T D Qi E^{6){E f + yi r 7 (l - *i) 1-7 S(*i,7). 



(B7) 



,K 



In the case of a flat surface the slant depth is expressed by R = h/ cos 9, where h is a vertical depth below the surface. 

It is important to note that formula( |B7[ ) evaluated in the assumption of the single slope (a, f3) of linear total losses 
a + (3E may be easily extended to the case of multislope losses. For example, in the case of a change of loss slope 
from (ai,/3i) to (0121P2) on the way of a muon along the path R the integral equation (2.3) is transformed to 



dE 



dE 



Ef a + PE J Ef ax+fixE J El2 a 2 + l3 2 E 



dE 



= R, 



(B8) 



where E\ 2 is energy value in the point of a slope change. Thus, to compute (2.2) one must determine E s from the 
integral equation 



dE 



E 12 OL 2 + (3 2 E 



= R-R 



12- 



where 



R12 



In 



Q-l + #12 & 

cti + Effii 



(B9) 



(BIO) 



It leads to the integral flux expressed by F c ;(> Ef,R,9) = F d {> E X2l (R — i?i 2 ),6*; a 2 ,P2)- In other words, when 
using expression (B7) for the slant depths R > R12 one must substitute R —* (R — R12) and Ef — » E12 and use the 
values (a 2 ,02) for a loss description. For slant depths R < R12 the use of (B7) remains unchangeable and the loss 
values are expressed by (ai,/?i). 

This algorithm may be extended to computations with any number of slopes of the energy dependence L(E). 
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FIG. 1. Effective cosine as a function of zenith angle. Curve (a) is geometrical effective cos ine cos#* given by Eq. (Al) of 
Appendix [a|. Curve (b) is effective cosine cos#** with the correction and is given by Eq. (3.4). Curves (c) and (d) represent 
the ratio cos#/ cos#* and cos#/cos#", correspondingly. 
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FIG. 2. Differential spectra of conventional muons at sea level for six zenith angles, sec 8: 1, 2, 3, 5, 10, and 20, from bo ttom 
to top. Curve labels correspond to values of sec 8. Solid curves - NSS spectrum, dotted curves - spectrum defined by Eq. (3.3), 
dashed curves - spectrum of Ref. [|16| with an angular distribution from Ref. km shown only for two values of sec 8 1 and 20. 
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FIG. 3. Sea level muon angular distribution defined by the ratio of differential flux at coaO to that at vertical direction. The 
disributions are shown for four values of energy Eo (TeV): 1, 2, 10, and 100, from bottom to top. Curve labels correspond to 
values of Eq. The regions of small (right picture) and big (lef t pi cture) zenith angles are enlarged for more details. Solid curves 



- NSS spectrum, dotted curves - spectrum defined by Eq. (3.3), dashed curves - are taken from Ref. [|25| (for tt + 15% K), 
and closed circles - spectrum of Ref. 11 51] . 
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FIG. 4. Ionization muon energy losses in pure water (as MUM code output). Dashed curve - a c given by classic Bethe-Bloch 
formula, solid curve - sum of a c and a e (where a e is loss due to bremsstrahlung with 7-quantum emission by atomic electron) . 
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FIG. 5. The b-terms of muon radiative energy losses in pure water (as MUM code output). Losses due to direct e + e~ pair 
production b p (curve (a)), bremsstrahlung bt (b), and photonuclear interaction b n (c) arc shown. Curve (d) is the sum of all 
b-terms. 
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FIG. 6. Muon energy losses in pure water as a function of energy E (as MUM code output). Curve (a) is loss due to 
ionization a{E) and (b) is total loss due to all radiative processes b(E)E. Curve (c) describes total energy losses a(E) + b(E)E. 
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FIG. 7. Survival probabilities in pure water P(Eo,R, > Ef ) as a function of starting muon energy Eo, resulted from MC 
simulations by using the MUM code as described in Sec, IV B (dots). Solid lines are drawn through the dots by the spline 
method. Four pictures are shown for various cut-off energies Ef. 10 GeV (a), 100 GeV (b), 1 TeV (c), and 10 TeV (d), 
correspondingly. Curve labels at each picture correspond to following set of slant depths R: 1.15 km (1), 2.07 km (2), 3.45 km 
(3), 5.75 km (4), 8.05 km (5), 10.35 km (6), 12.65 km (7), 14.95 km (8), 17.25 km (9), 19.55 km (10), 21.39 km (11), and 23.23 
km (12), from left to right. 
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FIG. 8. Survival probabilities P{Eo, R, > Ef) in pure water as a function of cut-off muon energy Ef for a given slant depth 
R — 1.61 km. The results of different MC codes of muon propagation are shown. Solid lines - MUM code, dashed - PROPMU 
code, circles - MUSIC code. The distributions are given for seven values of starting energy Eo (TeV): 0.5, 0.6, 0.7, 1, 3, 10, 
and 30, from left bottom to right top. 
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FIG. 9. The values of starting muon energy Eo at a sea level which result in passing till underwater vertical depth of 1.15 km 
at various zenith angles and cut-off energies Ef. Results obtained with MUM code with taking into account loss fluctuations 
are given. Solid curves are shown for Ef of 10 GeV (left picture) and 1 TeV (right picture). Dashed curves are shown for Ef 
of 100 GeV (left picture) and 10 TeV (right picture). Dotted curves are shown for Ef of 316 GeV at left picture. For each Ef 
we present the set of three curves corresponding to following meaning of Eq: E™ m , Eq and Eq 1<1x from bottom to top. E™™ is 
a minimal detectable starting energy Eo given by survival probability of 10 -4 . Eo is a mean sea level energy calculated using 
NSS spectrum. Eq 1<1x is that value of upper limit of the integral (1.5), which ensures 99.5% of underwater flux, calculated 
using NSS sea level spectrum with £™ aa; =36.3 PeV. 
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FIG. 10. The ratio of underwater integral flux expected at vertical depth /i=1.15 km in water at various zenith angles 
to that expected for vertical direction for the corresponding slant depths R — h/ cos9. Solid curves result from numerical 
computations allowing for loss fluctuations by using of NSS sea level spectrum. They are given for three values of cut-off energy 
Ef. 10 GeV, 1 TeV and 10 TeV, from bottom to top. Dashed curve is derived from sea level spectrum of Ref. [[l6| with an 
angular distribution from Ref. J2a] and is shown for _E/=10 GeV. 
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FIG. 11. Correction factor Cj as a function of slant depth R in pure water. The results obtained using sea level spectrum 
defined by Eq. (3.J) are given. Solid curves correspond to numerical calculations for vertical case 9 — 0°. Dashed curves 
describe the correction factor computed at vertical depth h of 1.15 km for various zenith angles as a function of slant depth 
defined by R = h/ cos 6. Both solid and dashed curves are shown for four values of cut-off energy Ef. 10 GeV, 100 GeV, 1 
TeV, and 10 TeV, from top to bottom. 
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FIG. 12. Correction factor C/ as a function of spectral index 7 of sea level spectrum for various depths in water for vertical 
direction. The distributions for cut-off ener gy Ef=10 GeV are given. Solid curves correspond to numerical computations by 



using sea level spectrum defined by Eq. (3.3) with varying spectral index 7. Open circles correspond to numerical computations 
using "simplified" VZK sea level spectrum, closed circles - Gaisser's sea level spectrum, squares - MACRO sea level spectrum 
(see Table |). All distibutions are shown for the following values of vertical depth in pure water: 1.15 km (a), 3 km (b), 5 km 
(c), 7 km (d), 9 km (e), 11 km (f), 13 km (g), 15 km (h), 17 km (i), and 21 km (j), from top to bottom. 
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FIG. 13. Underwater integral muon flux allowing for loss fluctuations as a function of zenith angle at different vertical 
depths. Four pictures are shown for various cut-off energies Ef. 10 GeV (a), 100 GeV (b), 1 TeV (c), and 10 TeV (d), 
correspondingly. Four curves at each picture correspond to vertical depths h: 1.15 km, 1.61 km, 2.0 km, and 3.0 km, from 
top to bottom. Solid cur ves result from numerical computations by using the NSS sea level spectrum and dotted ones from 
analytical expression (5.3) by using the sea level spectrum (3.3). 
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FIG. 14. Underwater vertical integral muon flux allowing for loss fluctuations as a function of cut-off energy Ef for various 
sea level spectra presented in Table rtl Two set of curves are given for two vertical depths R: 10 km (a) and 15 km (b). Dashed 
curves result from numerical computations by using of VZK sea level spectrum jlq], dotted - Gaisser's sea level spectrum, 
solid - par ame trization (3.2) of present work, dash-dotted - LVD sea level spectrum. Closed circles result from analytical 
expression (5.3) with using the same sea level spectra. 
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Underwater differential muon spectrum allowing for loss fluctuations as a function of energy Ef at vertical depth 
The curves are given for twelve zenith angles cos#: 1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 02, 0.15, 0.1, and 0.05, 



FIG. 15. 
of 1.15 km. 

from top to bottom. Solid curves result from numerical differentiation of integral flux expressed by (4.5) and based on NSS 
sea level s pect rum. Dotted curves result from numerical differentiation of the analytical expression (5.3) based on the sea level 
spectrum (3.3). 
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